function [DivorceCostsM,DivorceCostsF,DivorceCostsMF] = get_divorce_costs(singlefemale,singlemale,couple,leisurem,leisuref,hworkm,hworkf,consumption,wagem,wagef,nonlabor,ub)
%%
% LP to get the divorce costs (stability indices)
% Program uses YALMIP (https://yalmip.github.io) connected to CPLEX or GUROBI

% INPUT VARIABLES:
% singlefemale[households]  -- indicator for a single female household
% singlemale[households]    -- indicator for a single male household
% couple[households]        -- indicator for a couple household
% leisurem[households]      -- leisure per week of the male
% leisuref[households]      -- leisure per week of the female
% hworkm[households]        -- time spend on home production by the male
% hworkf[households]        -- time spend on home production by the female
% consumption[households]   -- weekly Hicksian consumption
% wagem[households]         -- wage rate of the male
% wagef[households]         -- wage rate of the female
% nonlabor[households]      -- weekly nonlabor income of the household
% ub                        -- upper bound on economies of scale of leisure (aleisure)

% OUTPUT VARIABLES:
% DivorceCostsM  -- Divorce Costs for individual rationality of males
% DivorceCostsF  -- Divorce Costs for individual rationality of females
% DivorceCostsMF -- Divorce Costs for no blocking pairs 

%%

%defining constants:
workh = 112;        % total potential working hours
num_of_couples = length(find(couple)); %number of couples in the data set
num_of_hh = length(couple); %total number of households in the data.

%% defining the decision variables

% private consumption 
privateqM = sdpvar(num_of_hh,1);                      % private consumption of Hicksian good by male
privateqF = sdpvar(num_of_hh,1);                      % private consumption of Hicksian good by female
privatemhworkm = sdpvar(num_of_hh,1);                 % private consumption of male's time spent on household work by male 
privatemhworkf = sdpvar(num_of_hh,1);                 % private consumption of male's time spent on household work by female
privatefhworkm = sdpvar(num_of_hh,1);                 % private consumption of female's time spent on household work by male
privatefhworkf = sdpvar(num_of_hh,1);                 % private consumption of female's time spent on household work by female
privateleisurem = sdpvar(num_of_hh,1);                % private consumption of male's time spent on leisure (only consumed by male)  
privateleisuref = sdpvar(num_of_hh,1);                % private consumption of female's time spent on leisure (only consumed by female)

% Barten scales for couples
aleisurem = sdpvar(1);                                % degree of publicness in leisure of male for couples
aleisuref = sdpvar(1);                                % degree of publicness in leisure of female for couples
ahworkm = sdpvar(1);                                  % degree of publicness in domestic work of male for couples
ahworkf = sdpvar(1);                                  % degree of publicness in domestic work of female for couples
aq = sdpvar(1);                                       % degree of publicness in Hicksian good for couples   

% Barten scales for singles
bleisurem = sdpvar(1);                                % degree of publicness in leisure of male for singles 
bleisuref = sdpvar(1);                                % degree of publicness in leisure of female for singles 
bhworkm = sdpvar(1);                                  % degree of publicness in domestic work of male for singles
bhworkf = sdpvar(1);                                  % degree of publicness in domestic work of female for singles
bq = sdpvar(1);                                       % degree of publicness in Hicksian good for singles

nl_incomeM = sdpvar(num_of_hh,1);                     % non-labor income of male
nl_incomeF = sdpvar(num_of_hh,1);                     % non-labor income of female

% stability indices
divorce_costsM = sdpvar(num_of_couples,1); % divorce costs are defined for every couple given every possible outside option including married and single potential partners
divorce_costsF = sdpvar(num_of_couples,1); % divorce costs are defined for every couple given every possible outside option including married and single potential partners
divorce_costsMF = sdpvar(num_of_hh,num_of_hh,'full'); % divorce costs are defined for every couple given every possible outside option including married and single potential partners


%% setting up the constraints
Constraints = [];

% Private consumption of Hicksian good must be balanced 
for i=1:num_of_hh
    if (couple(i) == 1)
    Constraints=[Constraints, ...
        consumption(i)*(1-aq) == privateqM(i) + privateqF(i), privateqM(i) >=0, privateqF(i) >= 0];
    else
    Constraints=[Constraints, ...
        consumption(i)*(1-bq) == privateqM(i) + privateqF(i), privateqM(i) >=0, privateqF(i) >= 0];
    end
end

for i=1:num_of_hh
      if singlefemale(i) == 1
            Constraints = [Constraints, ...
                privateqM(i) == 0];
        else
            if singlemale(i) == 1
               Constraints = [Constraints,...
                  privateqF(i) == 0];
            end
      end 
end

% Domestic work by male 
for i=1:num_of_hh
    if (couple(i)==1)
    Constraints=[Constraints, ...
        hworkm(i)*(1-ahworkm) == privatemhworkm(i) + privatemhworkf(i), privatemhworkm(i) >= 0, privatemhworkf(i) >= 0];
    else
        Constraints=[Constraints, ...
        hworkm(i)*(1-bhworkm) == privatemhworkm(i) + privatemhworkf(i), privatemhworkm(i) >= 0, privatemhworkf(i) >= 0];
    end
end

for i=1:num_of_hh
      if singlefemale(i) == 1
            Constraints = [Constraints, ...
                privatemhworkm(i) == 0];             % since hworkm is 0 for single female, this means both privatemhworkm and privatemhworkf are 0 in this case 
      else
            if singlemale(i) == 1
               Constraints = [Constraints,...
                  privatemhworkf(i) == 0];
            end
      end 
end

% Domestic work by female
for i=1:num_of_hh
    if (couple(i)==1)
    Constraints=[Constraints, ...
        hworkf(i)*(1-ahworkf) == privatefhworkm(i) + privatefhworkf(i), privatefhworkm(i) >=0 , privatefhworkf(i) >= 0];
    else
    Constraints=[Constraints, ...
        hworkf(i)*(1-bhworkf) == privatefhworkm(i) + privatefhworkf(i), privatefhworkm(i) >=0 , privatefhworkf(i) >= 0];
    end
end

for i=1:num_of_hh
      if singlefemale(i) == 1
            Constraints = [Constraints, ...
                privatefhworkm(i) == 0];
        else
            if singlemale(i) == 1
               Constraints = [Constraints,...
                  privatefhworkf(i) == 0];        % since hworkf is 0 for single male, this means both privatefhworkm and privatefhworkf are 0 in this case 
            end
      end 
end

% Leisure male
for i=1:num_of_hh
    if(couple(i)==1)
    Constraints=[Constraints, ...
        leisurem(i)*(1-aleisurem) == privateleisurem(i)];
    else
   Constraints=[Constraints, ...
        leisurem(i)*(1-bleisurem) == privateleisurem(i)];
    end
end

% Leisure female
for i=1:num_of_hh
    if (couple(i)==1)
    Constraints=[Constraints, ...
        leisuref(i)*(1-aleisuref) == privateleisuref(i)];
    else
    Constraints=[Constraints, ...
        leisuref(i)*(1-bleisuref) == privateleisuref(i)];
   
    end
end

% balanced nonlabor income:
for i=1:num_of_hh
        if couple(i) == 1
            Constraints = [Constraints, nonlabor(i) == nl_incomeM(i)+nl_incomeF(i)];
        else
            if singlefemale(i) == 1
                Constraints = [Constraints, nonlabor(i) == nl_incomeF(i)];
            else
                if singlemale(i) == 1
                Constraints = [Constraints, nonlabor(i) == nl_incomeM(i)];
                end
            end 
        end
end

% bounding the possible sharing of nonlabor income between partners [.4,.6]
for i=1:num_of_hh
    % check that dealing with couple
    if couple(i) == 1
        if nonlabor(i)>=0
            Constraints = [Constraints,
                .4*nonlabor(i) <= nl_incomeM(i) <= .6*nonlabor(i), .4*nonlabor(i) <= nl_incomeF(i) <= .6*nonlabor(i)];
        else 
            Constraints = [Constraints,
                .6*nonlabor(i) <= nl_incomeM(i) <= .4*nonlabor(i), .6*nonlabor(i) <= nl_incomeF(i) <= .4*nonlabor(i)];
        end
    end
end

% INDIVIDUAL RATIONALITY:
for i=1:num_of_hh
        if (couple(i)==1)
        Constraints = [Constraints,
            (wagem(i)*workh)*divorce_costsM(i) + nl_incomeM(i) <= privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(i) + privateqM(i) + hworkm(i)*ahworkm*wagem(i) + hworkf(i)*ahworkf*wagef(i) + aq*consumption(i) + leisurem(i)*aleisurem*wagem(i) + leisuref(i)*aleisuref*wagef(i)];
        % IR-W
        Constraints = [Constraints,
            (wagef(i)*workh)*divorce_costsF(i) + nl_incomeF(i) <= privateleisuref(i)*wagef(i) + privatemhworkf(i)*wagem(i) + privatefhworkf(i)*wagef(i) + privateqF(i) + hworkm(i)*ahworkm*wagem(i) + hworkf(i)*ahworkf*wagef(i) + aq*consumption(i) + leisurem(i)*aleisurem*wagem(i) + leisuref(i)*aleisuref*wagef(i)];
        end
end

% NO BLOCKING PAIR CONDITIONS
for i=1:num_of_hh
    for j=1:num_of_hh
        if (couple(i)== 1) && (couple(j)==1)
            Constraints = [Constraints, 
                (wagem(i) + wagef(j))*workh*divorce_costsMF(i,j) + nl_incomeM(i) + nl_incomeF(j) <= privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(j) + privateqM(i) + privateleisuref(j)*wagef(j) + privatemhworkf(j)*wagem(i) + privatefhworkf(j)*wagef(j) + privateqF(j) + max(hworkm(i),hworkm(j))*ahworkm*wagem(i) + max(hworkf(i),hworkf(j))*ahworkf*wagef(j) + aq*max(consumption(i),consumption(j)) + max(leisurem(i),leisurem(j))*aleisurem*wagem(i) + max(leisuref(i),leisuref(j))*aleisuref*wagef(j)];
        end
        if (couple(i)==1) && (singlefemale(j)==1)
           Constraints = [Constraints, 
           (wagem(i) + wagef(j))*workh*divorce_costsMF(i,j) + nl_incomeM(i) + nl_incomeF(j) <= privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(j) + privateqM(i) + privateleisuref(j)*wagef(j) + privatemhworkf(j)*wagem(i) + privatefhworkf(j)*wagef(j) + privateqF(j) + hworkm(i)*ahworkm*wagem(i) + hworkm(j)*bhworkm*wagem(i) + hworkf(i)*ahworkf*wagef(j)+ hworkf(j)*bhworkf*wagef(j) + aq*consumption(i) + bq*consumption(j) + leisurem(i)*aleisurem*wagem(i) + leisurem(j)*bleisurem*wagem(i) + leisuref(i)*aleisuref*wagef(j)+ leisuref(j)*bleisuref*wagef(j)];
        end
       if (couple(i)==1) && (singlemale(j)==1)
                    Constraints = [Constraints, 
           (wagem(j) + wagef(i))*workh*divorce_costsMF(j,i) + nl_incomeM(j) + nl_incomeF(i) <= privateleisurem(j)*wagem(j) + privatemhworkm(j)*wagem(j) + privatefhworkm(j)*wagef(i) + privateqM(j) + privateleisuref(i)*wagef(i) + privatemhworkf(i)*wagem(j) + privatefhworkf(i)*wagef(i) + privateqF(i) + hworkm(j)*bhworkm*wagem(j) + hworkm(i)*ahworkm*wagem(j) + hworkf(j)*bhworkf*wagef(i)+ hworkf(i)*ahworkf*wagef(i) + bq*consumption(j) + aq*consumption(i) + leisurem(j)*bleisurem*wagem(j) + leisurem(i)*aleisurem*wagem(j) + leisuref(j)*bleisuref*wagef(i) + leisuref(i)*aleisuref*wagef(i)];
        end
    end
end
   

% bounding the divorce costs
% individual rationality
for i=1:num_of_hh
    if (couple(i)==1)
    Constraints = [Constraints, 0 <= divorce_costsM(i) <= 1, 0 <= divorce_costsF(i) <= 1];
    end
end

% no blocking pairs
for i=1:num_of_hh
    for j=1:num_of_hh
         Constraints = [Constraints,  0 <= divorce_costsMF(i,j) <= 1];
    end
end

for i=1:num_of_hh
    for j=1:num_of_hh
        % bounds on divorce costs work only for m_i and w_j who are in the consideration set of each other
         if (couple(i) == 1) && (singlemale(j) == 1)
             Constraints = [Constraints, divorce_costsMF(i,j) == 0];
         end
         if (singlefemale(i) == 1) && (couple(j) == 1) 
             Constraints = [Constraints, divorce_costsMF(i,j) == 0];
         end
         if (couple(i) == 0) && (couple(j) == 0)
             Constraints = [Constraints, divorce_costsMF(i,j) == 0];
         end
         Constraints = [Constraints,  0 <= divorce_costsMF(i,j) <= 1];
         
    end
end

Constraints = [Constraints,  0 <= aq <= 1, 0 <= aleisurem <= ub, 0 <= aleisuref <= ub, 0 <= ahworkm <= 1, 0 <= ahworkf <= 1];
Constraints = [Constraints,  0 <= bq <= 1, 0 <= bleisurem <= ub, 0 <= bleisuref <= ub, 0 <= bhworkm <= 1, 0 <= bhworkf <= 1];


%% Specifying Objective Function
Objective = sum(divorce_costsM)+ sum(divorce_costsF)+ sum(divorce_costsMF(:));

% Solving the problem
options = sdpsettings('verbose',0,'solver','gurobi');

% Solve the problem
% -Objective because the default is to minimize while we want to maximize.
sol = optimize(Constraints,-Objective,options);

% Analyze error flags
if sol.problem ~= 0
 display('Hmm, something went wrong!');
end

%getting the value out
DivorceCostsM  = value(divorce_costsM);
DivorceCostsF  = value(divorce_costsF);
DivorceCostsMF = value(divorce_costsMF);

return
